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<N ■ Abstract 

N 

The present review is devoted to our recent studies on the excitonic motion in photosynthetic 
systems. In photosynthesis, the light photon is absorbed to create an exciton in the antenna 
complex of the photosynthetic pigments. This exciton then migrates along the chain-biomolecules, 
O ■ like FMO complex, to the reaction centre where it initiates the chemical reactions leading to 

biomass generation. Recently, it has been experimentally observed that the exciton motion is 

highly quantum mechanical in nature i.e., it involve long time (~ 600 femto sec) quantum coherence 
-i— > ' 
cZ2 ' 

effects. Traditional semiclassical theories like Forrester's and second Born master equations cannot 

be applied. We point out why the 2nd Born non-Markovian master equation and its Markovian limit 



*^3 ' (also called the Redfield master equation) cannot be used to explain the observed long coherences. 

Briefly, the reason is that these approaches are perturbative in nature and in real light harvesting 



systems various couplings (system-system and system-bath) are of the similar order of magnitude. 



£> \ Various new approaches are being developed to go beyond the above two limiting theories. The 

re- 
present review is not a review in the usual sense of the word as we summarize our own approaches 

and only refer to the literature for the other ones. A brief introduction to the sophisticated 2D 

^i . photon echo spectroscopy is also given at the end with an emphasis on the underlying physics of 



the multidimensional echo spectroscopies. 
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BChl molecules 



FIG. 1: On the average, on earth, biomass worth twice the mass of the Great Pyramid of Giza 
(~ 10 million tons) is being produced every hour. (Right) FMO protein acts as a "wire" which 
connects an antenna to a reaction center. 

I. INTRODUCTION TO EXCITATION ENERGY TRANSFER (EET) AND THE 
STATEMENT OF THE PROBLEM 



Photosynthesis provides chemical energy for almost all life on earth. Understanding of 
the natural photosynthesis can enable us to construct artificial photosynthesis devices and 
a solution to the future energy problems. We know that the safe (green) energy resources 
is a big challenge of the future and we know that our present energy technology is seriously 
disturbing our environment [1]. This prompts us to investigate " Green" resources of energy 
and photosynthesis is one of them. Photosynthesis is an interesting phenomenon. On the 
average, on earth, biomass worth twice the mass of the Great Pyramid of Giza (~ 10 million 
tons) is being produced every hour or ~ 3 million kg/sec. 

In photosynthetic systems a central role is played by the energy transport "wire" -the 
FMO protein which is a trimer made of identical subunits containing seven bacteriochloro- 
phyll (BChl) molecules each (Fig. 1). The photosynthesis process can be divided into the 
following steps [2]: 

1. A light photon is absorbed to create an excitation in the antenna complex of the 
photosynthetic pigments. 

2. This excitation then migrates along the chain-biomolecules, like FMO complex, to the 
reaction center where it initiates the chemical reactions leading to biomass generation. 

Quantum dynamics of EET (excitation energy transfer) can be very easily analyzed in 
the following two limiting cases. We identify first the couplings: 
Two important couplings: 



Incoherent Case 



Coherent Case 



FIG. 2: Two extreme cases: (1) Upper: incoherent (classical diffusive motion), (2) Lower: ultra- 
quantum (delocalized excitation). The theoretical investigation of the intermediate case is a big 
challenge. 

1. Inter BChl molecules (system-system) coupling J (which is responsible for EET). 

2. The system-bath (BChl molecule-protein) coupling A (which is responsible for deco- 
herence) . 



Two important timescales: 

1. Excitation transfer time scale T trans f er = 4 ~ 265/s, for J ~ 20cm _1 ~ 4 x 
10~ 22 joules (usual in photosynthesis pigments). 

2. Decoherence time scale Td eC o — f • 

If the system-bath coupling is very weak and r^eco » Ttransfer, the system is almost 
closed and dynamics is quantum mechanical in nature, i.e., one can use the Shroedinger's 
equation to analyze it in the extreme case. But in the opposite case r^ eco « Transfer 
(strong system-bath coupling), the system is almost open, decoherence rate is very fast and 
the dynamics is almost incoherent. One can analyze the process with simple Pauli type 
master equation with rate of transfer of the excitation from one molecule to other calculated 
with Fermi's golden rule (Forester's theory) [3]. 

The important problem arises in the intermediate regime as real light harvesting systems 
do not fall in either of the extreme cases. Recently, it has been experimentally observed that 
the exciton motion is highly quantum mechanical in nature i.e., it involve long time (~ 600 



fs) quantum coherence effects[4]. These discoveries caused a lot of excitement in field[5] as 
the traditional view was of incoherent "hopping" of the exciton(Fig. 2). In the weak system- 
bath coupling case, the standard approach was the 2nd Born quantum master equation which 
is a perturbative quantum master equation (upto second order in system-bath interaction). 
Its Markovian and secular approximation is known as Redfield master equation[6]. The 
2nd Born quantum master equation can be obtained from the Nakajima-Zwanzig projection 
operator technique by restricting the perturbation series upto second order [6]. These simple 
and powerful projection operator techniques were introduced by Zawnzig in 1960's in then 
active field of non-equilibrium statistical mechanics. As the the light harvesting pigments fall 
in the intermediate regime (system-bath coupling is of the order of system-system coupling) 
one clearly cannot use the 2nd order perturbative quantum master equation for its study in 
its original form. But in its modified form its scope becomes wider [7]. 

It is of interest to quantitatively know upto what value of system-bath coupling strength 
and other important couplings in the problem, one can use the 2nd Born master equation. 
Section II deals with this. 

In an important case of fast bath relaxation (when bath degrees-of-freedom re-organize 
very fast as compared to the transfer time scale of the exciton) a very useful approxima- 
tion can be made. The details of which are given below. This is called the Markovian 
approximation. In the following sections we summarize our study [8] of the quantitative de- 
termination of the regime of validity of the second order approximation and the Markovian 
approximation. 

II. MICROSCOPIC APPROACH: 2ND BORN MASTER EQUATION 

As is well known that the 2nd Born quantum master equation can be obtained from the 
Nakajima-Zwanzig projection operator technique by restricting the perturbation series upto 
second order in the system-bath interaction [6]. In the following we will apply this master 
equation to a concrete model of a dimer (open two-state quantum system) which caricature 
the dynamics of decoherence in a typical photosynthetic system[9]. 



A. Non-Markovian solution 

Projection super-operators and dynamics of the relevant system: The total system (rel- 
evant system (electronic part) + Bath (phonons)) dynamics is pure quantum in nature. 
The partial time derivative of the total density matrix is given by Liouville-van Neumann 
equation (classical equivalent is the invariance of the "extension" in phase space): 

^%^ = C(t)fLa(t) = -{[HLvu + Hi plat)] (1) 

Here the interaction representation is used J (t) = U s (t)OUs(t), Us(t) = Exp(—j^H s t) 
(see for details any standard refs[6]). H e i_ p h is the system-bath interaction Hamiltonian 
(for the 2nd Born approximation ||.H e |_p/ l || <C ||-ff s ||- H\, is the bath Hamiltonian. To 
construct the equation-of-motion for the relevant part of p[ ota iit) i.e., p(t), one defines the 
super-operator: 

VO = R eq tr R (6), (2) 

called the projection super-operator, one also defines Q = I — V. By applying V, Q on 
the Liouville-van Neumann equation turn- by-turn, we get: 

dVpIt ^ l{t) =VC{t){V+Q)p\ otal {t) 

dQ ^ al{t) = Qm(V+Q)pLai(t). (3) 

Solving the second equation formally for the irrelevant part (QPtotaii^)) > an d inserting 
in the first, one obtains the required equation-of-motion for the relevant part (p(t)R eq = 

VphJf)) ( see for details [6]): 

With 2nd Born approximation (i.e., by expanding the time evolution operator upto the 
first power in the system-bath interact ion [6]) and for a traditional dimer system [9]: 



H tot = H s + H ph + H el -P h 

H s = H el + H reor9 

2 

H el = 5»>H + J(|1><2| + |2}(1|) 

n=l 

2 

n=l i 

2 

n=l i 

2 

H el ~ ph = J^K^n, K = |n)(n|, u n = -^2hwid ni q h (4) 

n=l i 

master equation takes the form, 



<V(£) i ^ 7 7 



1 E f^t-rj^W.WW] -C^(t-r)^(0y(r)^/(r)]) (5) 

Here in the Hamiltonian, |n) represents the state in which ONLY nth site is excited and 
all others are in the ground state i.e., \n) = |0 n ,e)|0m^n, 5 )- H s is the system Hamiltonian 
which consists of H el the electronic Hamiltonian for the two level system, and H reor9 the 
Hamiltonian for the re-organization energy (the elastic energy related to the physical orga- 
nization of the bath degrees-of- freedom) . H ph is the phonon Hamiltonian and H el ~ ph is the 
system-bath coupling Hamiltonian. In the absence of phonons, e° is the excited electronic 
energy of n th site and J is the electronic coupling between the sites which is responsible 
for excitation transfer. The ground state energies of both donor and acceptor are set equal 
to zero and Xj is the re-organization energy of the j th site (Dissipated energy in the bath 
after the electronic transition occurs). dji,qi,Pi are the dimensionless displacement of the 
equilibrium configuration of the i th phonon mode, dimensionless coordinates, momenta of 
the i th phonon mode respectively. 

In the master equation, the bath correlation functions (bath is assumed to be a continuum 
of harmonic oscillators (valid when an-harmonic terms are not important)) are: 



<?«(*) = <«#)«,(())) -<«*)<«*) (6) 

We consider a case where the characteristics of the bath as seen by both the sites are 
the same, and there is no systematic bath correlations between the sites. Thus the bath 
correlation function takes the form: C{j(t) = C{t)8if 

C(t) = I ^C(u)e-^. (7) 



oo 



C(w)=2h(l + n(u))J(w), J(u) = 2A ^ _ (8) 

ur + 7^ 

Assuming the Drude-Lorentz model[9] for the bath spectral density, and assuming the 
high temperature approximation {-rHp « 1), as appropriate for the FMO problem, we 
obtain 

o\ 1 

c ^-r *• "=fer < 9 > 

B. Representations 

The equation (5) is an operator equation and this can be expressed in site or in energy 
representation. In site representation, with definitions x(t) = pn(t) = (1|/3(£)|1) (site), 
7/i (t) = Re[pi 2 (t)], and y 2 (t) — I m [pi2(£)], an d with lengthy but straightforward calculations 
(see for details [8]), the equation (5) can be written explicitly as a set of coupled integro- 
differential delay equations: 

ij f- ->> 

[771 cos(E 12 (t - r))j/i(r) + r/ 2 sin(E 12 (t - r))y 2 (r)] 

[-r/ 2 sin(E 12 (t - r))7/i(r) + 7/3 cos(£i 2 (£ - r))7/ 2 (r) + 2fij/ 2 (r)] 
with t/! = 1, r/2 = - 7S A_, 7/3 = ^,, £ 12 = (^ - £ 2 )/ft = -^S « = 

2J2 

A 2 +4J 2 • 



For energy representation we need the eigensystem of the Hamiltonian. Let the kets 
|ei,2) be the eigenstates of the Hamiltonian H s . The reduced density matrix in energy 
representation can be expressed as: 

Pl b (t) = (e a \p(t)h), (10) 

with time evolution given as, 
dp e ab (t) 



1 f l 

? E J ^(i 



-tUabPab - T? 7 , / ar( ^ { t ~ T 



i,c,d=l J ° 



dt ao,ab h 2 

[yacycd e -iu, cb (t-T) p eJ^ _ yacydb e ^ ad {t-r) p e cd ^ 

-C*(t - T) lVrV? b e-^t-r) p e cAT) _ ycdyd^-i^it-r)^^ ? (n) 

with uj a ij = (E a — E b )/h and 

yac _ ^a^c y ac _ (12) 

1 v^TTv^f+T' 2 v^f+iy^+T' 

Eigensystem of the Hamiltonian: Assuming Ai = A2 = A the eigenvalues Ei and eigen- 
vectors \ej) of the system Hamiltonian 

2 



# s = ;O e n + A *>H + </(ii><2i + i2><ii) 



71=1 

can be easily obtained as 
1 



Ei,2 = ^(e? + e° + 2A =f y/(4 + e° + 2\f - 4(e?e° - J 2 + A(e? + e°) + A 2 )) 

l e l) = / o = 1 ' l e 2/ 



v^TTViy v^fTTVi 



ai, 2 = — (A=fVA 2 + 4J 2 ), A = e°-e°. 

Here the column vectors denote components in the site basis. The eigenkets are normalized 
and are orthogonal ai«2 = — 1- 

C. Numerical Approach 

It is well known that the numerical propagation of integro-differential equations is an 
involved task and time consuming. In the following we construct a simple method to the 
solution of numerical integration [8] . Specifically we utilize the exponential nature of the bath 



correlation function which helps to convert the set of coupled integro-differential equations 
to a bigger set of ordinary differential equations: 



h(t')= I dr'e' 

Jo 



cos[— (t' - r')]yi(r) + 7/2 sin[— (£' - r')]y 2 (T') 

7 7 



h{t') 

/s(0 



e T ' 'y 2 (T')dr' , 



dr'e 7 



77> Tp 

-77 2 sm[— (t' - r , )]3/i(r') + r7 3 cos[— (£' - r')]^2(r') 

7 7 



(13) 



With £' = 7t, t' = ■yr'. Here, t/i(t') = yi(t'/j), y2(t') = 2/2(^/7) and we also define 
x{1)=x{t/i). 

We obtain a set of coupled ordinary differential equations (note that these are much 
simpler to solve as compared to coupled integro-differential equations): 

1J 

x{t') = -y 2 (t'), 

*yn 






(3-f 2 H 2 



ne-*'Mt>)--^e-t'f 3 (t'), 



(3j 2 h 2 



£, 



E 



12 



7 



A(0, 



fi(t') - e*Ut) = e*M1) + — e'^Ut') - 

7 

Mt') = e«ifc(t), 



(14) 



7 V 7 

Our aim is to use this to establish the parameter range over which the Markovian approx- 
imation is valid. Before doing so we compare this method with the traditional method of so- 
lution of integro-differential equations[10]. In the straightforward numerical method (tradi- 
tional method) the integro-differential equation is first written as integral equation with dou- 
ble integration as dx{t)/dt = f Q f(x(t—r),t)dr converted to x(t) = j Q di! f Q drf(x(t'—T),t'). 
The double integration is then done self-consistently with numerical integration[10]. 

Speed: On the Lenovo ThinkCentre-i7, the traditional method took ~ 15 minutes to 
obtain the result, whereas the present method took only about a few milliseconds. A sample 
comparison is given in Fig. (3). 
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FIG. 3: Comparison of numerical results for a traditional method (blue-solid curve), to that in- 



troduced here (red-dotted curve). Parameters used are: 7 
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"\ A = 2 cm" 1 , J 



50 cm , A = 100 cm . Abscissa is in fs (fs = femto seconds). 



D. Markovian limit 

As mentioned before in an important case of fast bath relaxation (when bath degrees- 
of-freedom re-organize very fast as compared to the transfer time scale of the exciton) a 
very useful approximation can be made. This is called the Markovian approximation [21]. 
To do this approximation, note that it is particularly simple to invoke it in the energy 
representation (as one can "average out" the fast oscillations in the system's density matrix 
and can compare the temporal envelope of the system's density matrix with time decay 
of the bath correlation function). Hence, below we first utilize the energy basis and then 
convert the result back to the site representation. 

The Markov approximation can be performed when the time scale on which the envelope 
of the density matrix decays is much longer than the decay time of the phonon correlation 
function [6] . One can then introduce the following approximation: 



P £ ab (t-r) 



-iu) ab (t-r) ~e 



Pl b (t ~ T) 



-iu ab (t-T) ~e 



~pl b (t) = e^p e ab {t)- 



To do this energy representation equations are first converted to dimensionless form 
with t' = 7T (see for details [8]). Putting t — r = r' in the resulting equations and then 
implementing the above approximation on the density matrix elements allows the time 
integration to be performed easily for the case of exponential phonon correlation function. 
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The result is the set of Markovian equations: 



Kbit') = -iuab~p e ab {t') 

i,c,d x 

2A ^-^ / V ac V db , V cd V- 



/ yacydb ycdydb \ 

r-J \ 1 - «w ca 1 + iuj cd ) 



,2 
I .(■ 



Here p^ b {t' /j) = pl b (t'), u a b = co a b/j- The results can then be transformed back to the site 
representation using p^t) = (z|p(t)|j) = Ea,b(*l e a)Paft( e &li)- 

E. Limiting Cases: Analytical Results 

Markovian and non-Markovian results were obtained computationally and compared for 
various regimes. Before presenting them we show some interesting analytic results in two 
extreme cases on the parameter dependence of the region of validity of the Markov approx- 
imation: 

1. Strong Coupling Case: J> A 

For J ^> A, we have a\ — 1, «2 — — 1, and Vl° ~ 1/2 for i = j and ~ —1/2 for 
i y^ j ({i,j} = 1, 2) and V£' 3 — 1/2 for all i,j. One can then analytically solve the coupled 
equations to obtain the simple expression 

Pli(t') = -(e-^W)' +1), (15) 

for the traditional initial conditions p\i(t' = 0) = 1, p^f^' = 0) = p\i(t' = 0) = 0. The 
Markov approximation can be performed when the time scale on which the envelope of the 
density matrix decays is much longer than the decay time of the phonon auto-correlation 

function. Hence, 



4A 

< 1 



/3(4J 2 + frV) 



in the J ^$> A domain. 



2. Weak Coupling Case: JcA 



must hold for the Markov approximation to be valid 



For this case, we have ot\^ = o7(A =F a/A 2 + 4J 2 ) ~ ^/(l =p 1). Hence, in this domain 



2JV" ^ v " ^" ) — 2J V 

n T/22 ^ i T/ll ~ 1 T, 



ai ~ 0, and a 2 ~ A/ J. This leads to V^ 11 = V ± 12 = V 21 ~ 0, Vf 2 ~ 1. V 2 n ~ 1, K, 12 
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V 21 ~ J/A, and V 22 = (J/ A) 2 . 



Pn(t') 

2/A 

H 2 P<y 2 



(2(j/A) 2 (r + r*) - 4(j/A) 2 (r + r*)^(t') + 2(j/A)(pf 2 (t') + &(*))) 



&(0 = ^(0 + 

+^ ((J/A)r(2^(0 - 1) - (i + 2r(j/A) 2 ) P i 2 (t') + 2(j/A) 2 r^(0) ■ 



where T = 1 a_ . See for details the second paper in [8]. By separating real and imaginary 

tvy 

parts as pl 2 (t') = x{t') + iy(t') and writing p^t') = r(f), we have: 

N - <Rcos(7/£')e" ? '' + [ar/ + fa^sin^i')^*'), 

x(t') = e" ?t '(acos(r/t / ) - 6sin(r/t')), 

y (t') = e"^'(asin(r/t') + bcos(rjt')). (16) 

with initial conditions r(£' = 0) = 1, x(t' = 0) = a, y(t' = 0) = b. Here 
£ = 4A/(/i 2 /37 2 ), 7? = A/ fry, and e = J/A. Thus, for the Markov approximation to hold 

These are summarized in Table I. A numerical verification of 



requires 



£ = 4A/ / 3fr 7 2 < 1 



these analytical inequalities is given in Fig. (4). 



TABLE I: Regimes of validity of the Markov approximation 
Case Markovian approximation 

J » A 0(4 A-y') « l 



J«A J^«l 



WJFj 



F. Computational Results 

To investigated the validity regime of Markovian approximation beyond the above two 
limiting cases we have to rely on the numerical computation (as the analytic approach 
becomes very cumbersome). We here display an extensive list of graphs (Figs. 5,6, and 7) 
which explore how Markovian approximation behaves for various values of parameters. 

First consider the non-Markovian results (solid curves). Coherent oscillatory dynamics up 
to substantial time scales are evident. Oscillatory behavior in the populations are seen (Fig. 
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FIG. 4: Sample verification of the analytic inequalities (Table I). Time evolution of population on 
site 1: blue (solid) curve is the non-Markovian solution and red (dotted) curve is the Markovian 
approximation. 
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FIG. 5: Time evolution of population on site 1 [pn(t)] and the coherences [pi2(£)] [blue (solid) 
curve is the non-Markovian solution and red (dotted) curve is the Markovian approximation], for 
various values of A (in cm ). Other parameter are: A = 100 cm" 1 , J = 50 cm -1 , 7 = 10 13 s -1 . 
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FIG. 6: Time evolution of population on site 1: blue (solid) curve is the non-Markovian and red 
(dotted) curve is the Markovian) for various values of the reorganization energy A and inter-site 
coupling J. 
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FIG. 7: Time evolution of population on site 1 [blue (solid) curve is the non-Markovian and red 
(dotted) curve is the Markovian] for various values of the reorganization energy A and phonon 
relaxation time 7 _1 . 

5) to be accompanied by oscillations in the off diagonal elements p\% representing coherences. 
The oscillations fall off faster with increasing re-organization energies, as expected. Another 
point to be noted is that the relaxation to equilibrium populations occurs on a longer time 
scale than does relaxation of the coherences to zero. This difference is more evident at 
smaller values of re-organization energy A. The dependence of population relaxation on the 
bath correlation decay time 7 _1 is shown in Fig. 7. Clearly, the larger the 7 (i.e., fast bath 
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FIG. 8: Approximate Physical Picture 

relaxation), the better is the Markov approximation. 

Figures 5-7 also contain a comparison of the Markovian limit to the non-Markovian solu- 
tion for domains other than those in Table I. For the standard electronic coupling parameter 
values in photosynthetic EET: 7" 1 = 100 fs, J = 50 cm" 1 , A = 100 cm" 1 , T = 300 K, 
figure 5 shows that the Markovian approx is very good for A = 1 cm -1 , fair for A = 2 cm -1 , 
and invalid for reorganization energies A > 10 cm -1 . As typical values of A in photosyn- 
thetic EET systems are considerably larger than A = 10 cm" 1 , these results support the 
conclusions of Ref . [9] , although with a different approach. 

Figs. 6 and 7 show the validity of the Markov approximation obtained by varying (J, A) 
and (7 -1 , A) respectively, but keeping A = 100 cm" 1 . The results show that the Markovian 
approximation is poor for large A and small J, and for large 7" 1 and small A. Other 
parameter values can be easily examined using this approach. 

From these qualitative conclusions a rough physical picture can be drawn as depicted 
in Fig. (8). In the Markovian regime, after photo-excitation, bath (nuclear co-ordinates) 
re-organize very fast representing an "apt" bath (right hand picture) while in the non- 
Markovian regime bath correlations exist for longer time scale (of the order of exciton transfer 
time scale). Thus the combined (system+bath) dynamics is much more complicated in the 
non-Markovian regime. 

Here, results are given for the particular initial conditions: pn(0) = x(0) = 1, ?/i(0) = 
2/2(0) = 0, /i(0) = ^2(0) = ^3(0) = /i(0)/ 3 (0) = 0. These initial conditions (corresponding 
to all the population being on site 1, and no coherences) are those which have been used 
extensively in literature[9] but are somewhat unphysical[ll], because they lack initial coher- 
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ences which become important in preparatory photo-excitation. We have considered this 
problem in the second reference of the list [8] which considers photo-excitation of a dimer 
oscillator system with an ultra-short laser pulse. It appeared that the presence of initial 
coherence ( at t' — 0) effected the time scale on which the populations reach equilibrium 
value but had little effect on the overall damping-out of the coherences (see for details [8]). 

III. PHENOMENOLOGICAL APPROACH: A STOCHASTIC MODEL 

In the previous sections we saw that 2nd Born quantum master equation cannot be 
applied to the real light harvesting systems as in these systems the system-bath coupling 
(A ~ 100 crrT 1 ) is of the same order of magnitude of the system-system coupling (J ~ 
100 cm~ l ). To make the situation more intractable, it is not possible to justify the Markovian 
approximation when A>1 crrT 1 (given that Markovian master equations are much easier 
to solve than the non-Markovian ones). Thus the use of Markovian Redfield theory to these 
systems is questionable as pointed out in [8, 9]. This open up a difficult problem. One should 
formulate some non-perturbative theories. Recently Ishizaki and Fleming [9] have developed 
a formalism which goes beyond the limitations of the 2nd Born master equation. They use 
the reduced hierarchy equation approach previously developed by Tanimura and Kubo[12]. 
There is an other route to the problem pioneered by people like Silbey[13]. In this approach 
one uses a unitary transformation (called polaron transformation) to completely eliminate 
the system-bath coupling Hamiltonian. But this re-normalize the system Hamiltonian. Then 
one re-partition the resulting system Hamiltonian to identify a weaker term which can be 
used as a perturbation. The remaining problem is done in line with 2nd Born master 
equation[14]. 

We have developed an alternative stochastic approach which is phenomenological in 
nature[15]. This approach, as its input, takes the homogeneous line width from the ex- 
periment and uses Kubo's stochastic theory of motional narrowing to get phenomenological 
decoherence rate. 
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FIG. 9: Two interacting molecules 
A. The model and its solution 

We again consider the dynamics of exciton transfer between two molecules modeled as 
two-level electronic systems (Fig. 9). These two-level systems are electronically coupled with 
each other with coupling J. 

Due to the electronic coupling between the molecules the exciton will transfer back- 
and-forth between the molecules. This will happen forever if the molecules are completely 
isolated — pure oscillatory quantum motion. Now consider that our two-molecular system is 
open i.e., interacting with the external bath degrees-of-freedom — -the phonons. It is well 
known that the dynamics remain quantum at short time scales and becomes classical at 
longer time scales[16]. This quantum-to-classical crossover happens at a critical time scale 
which is inversely proportional to system-bath coupling energy (t c oc -). Large r\ (system- 
bath coupling) means fast quantum-to-classical crossover and vice versa. We denote system- 
bath coupling strength with r\ in the subsequent subsections (before we used A). 

In a recent contribution[15] we extracted r\ from experimental information using Kubo's 
stochastic theory of line shapes and observed upto what timescale one could see quantum 
effects. 

The important point is that we modeled the dynamical effect of phonons as a stochastic 
noise. The total Hamiltonian takes the form 

H = ei |l)(l| + ( Ca + e(t))|2)(2| + J(|l)(2| + |2)(1|). (17) 

Here ei is the energy of the upper electronic level of the first molecule and €2 is that of the 
second molecule. The ground state energies of both the molecules are taken to be zero. The 
energy separation €2 — e\ has a random component (due to phonons) which we denote with 
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e(t). e(t) is a stochastic process taken here as Gaussian White Noise (GWN): 

(e(t)) = 0, < C (t)e(T)) = ftV(t-r). (18) 

Here, as mentioned before, rj is the strength of system-bath coupling (also known as dynam- 
ical disorder) measured in the units of frequency. 

We start with Liouville-von-Neumann equation for the total density matrix, 

ih^ = [H,p(t)}. (19) 

As H is a stochastic operator, thus p(t) is also a stochastic operator. Therefore, we need 
to evaluate the averaged density matrix. So we need to do an averaging over the dynamical 
disorder which is denoted by (...). We define <;(£) = (p(t)). 

Averaging over dynamical disorder: 

^^ = (£p(^-(pW£^- (20) 

TermI Termll 

In term I and II above we have terms like (p(t)e(t)). As p(t) is a functional of e(t) (a 
stochastic quantity) the p(t) will also be a stochastic function. To decouple these we will 
use the famous theorem of Novikov[17]: 

«t) Pab (t)) = jT dt'(e(t)e(t')) (j^) (21) 



$Pab(t) 

6e(t>) 

some simplification(see appendix A), we get 



Here -§fnjy is the functional derivative. Using the properties of stochastic noise and with 



(e(t)p l2 (t)) = ihr]q 12 {t). (22) 

Finally, one has a set of coupled differential equations: 
cfcu(£) 



dt 

(ki2(t) 



-i(J/h)(c; 21 (t)-^ 2 (t)) 

-i(A/h)q 12 (t) - i{J/H)^ 22 {t) - q n (t)) - ^ 12 (t) 



dt 
SLi(t)+fe2(*) = l. (23) 

The above system of ODEs can be solved analytically, however, the exact expression is very 
cumbersome. We give the analytic solution only in the long time limit (see Appendix B). 
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FIG. 10: (Left) Schematic line broadening information in 2-D photon echo spectrum shown without 
cross peaks. The linewidth due to homogeneous and in-homogeneous broadening are in orthogonal 
directions as shown. (Right) The energy levels of the model system 

To simulate the dynamics of decoherence in this dimer model with physical parameters 
of the FMO problem (A = e\ — e 2 ~ 100cm -1 , J ~ 100cm" 1 ), we need to find out our 
phenomenological parameter n, for this we use Kubo's stochastic theory of lineshapes. 

B. Kubo's stochastic theory of lineshapes and estimation of 77 from motionally 
narrowed lineshape 

We now determine the phenomenological parameter n. Let us focus on exciton 1 in 2-D 
photon echo spectra which occurs at 810 nm (Fig. 10). 

We use Kubo's randomly modulated oscillator model for the exciton under question. The 
energy levels of the model system are given in Fig. 10. The levels 1 and 2 are separated on 
the average by uq = 810 nm. 

Let the level 2 be randomly modulated with a random process ooi(t) such that 
liniT^oo ^ f Q u>i(t)dt = 0. The random frequency of the level 2 is then given as u(t) = 
Uq + 0Ji(t). We use the classical oscillator model for the molecule with resonance frequency 
u = 810 nm. The equation of motion is x{t) = iu(t)x(t), with solution 

x(t) = x(0) exp I icu t + i I u x (t')dt' J . (24) 

The absorption spectrum measures a large number of possible realizations of the random 
process cui(t) ( The radiation field acts on a macroscopic number of molecules present with 
different realizations of cui (t) in the sample) . Thus one need to consider an ensemble average 
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to compare with absorption spectrum: 

(x(t)) = x(0)e iaj ° 4 /exp (i f Ul (lf)dA \ , (25) 

and the average time correlation is given by (x(t)x*(0)) = \x(0)\ 2 e luJot <j)(t) where 

<f>(t) = /exp (i J cuxit^dt'] \ (26) 

is called the relaxation function of the oscillator. The absorption spectrum is 



1 /— t-oc 

I(u -u ) = — e-^-^^(t)dt. (27) 

27T J-oo 

This is the direct consequence of the famous fluctuation-dissipation theorem. The temporal 
character of the decay of fluctuations tells directly the dissipative characteristics of the 
system. The intensity distribution I(oo) will be broadened by the stochastic process coi(t). 
To define the stochastic process Ui(t) let P(ui)dui be the probability to find the random 
frequency 0J\ to be in the range U\ to oj\ + dui when picked randomly from the ensemble. 
In Kubo's theory, the stochastic process is defined with two parameters (1) magnitude of 
modulation A 2 = j ufP(ui)dui = (uf) and (2) correlation time r c = L f c (t)dt where 
correlation function is defined as / c (r) = -gz (ui(t)u)i(t + r)). 

It is well known that for a Gaussian process, relaxation function can be written in terms 
of the correlation function: 

</>(£) = exp (-A 2 J (t - r)f c {r)dr\ (28) 

In the present case we have considered Gaussian White Noise (GWN) which has zero cor- 
relation time. Our case corresponds to the fast modulation case of Kubo r c << ^[18]. 
The correlation function decays very fast and the upper limit of the integral in the above 
equation can be extended to oo. This leads to <p(t) oc e _A Tc '*L This results in the famous 
narrowing of the lineshape from the Gaussian to Lorentzian form. In the present case we 
have ( /^ /~) = V^( T ) which leads to <p(t) = e _r? '*'. Comparison with the previous (p(t) 
shows that rj = A 2 r c which gives the rate of decay of the correlation function. In our GWN 
case A — > oo since white noise contains all frequencies and r c — > (delta correlated noise) 
but A 2 r c is finite and is equal to rj. 

Thus the absorption spectrum takes the form 
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I(u-Uo) = -- f \— 2 - (29) 

This is the famous Lorentzian lineshape narrowed from the Gaussian shape (called mo- 
tional narrowing). 

Our aim is to find out our phenomeno logical parameter rj. Thus we need to fit this 
I(u — Uq) with the real experimental observation and to extract rj. We will use this to 
simulate the quantum dynamics of the density matrix elements. 

The basic problem with linear absorption line shape is that it is broadened both by 
homogeneous and in-homogeneous mechanisms. In our case the broadening is homogeneous 
due to dynamical disorder and thus we need to subtract the in-homogeneous component due 
to static disorder. But thanks to the 2-D photon echo spectroscopy one has the important 
information about both homogeneous and in- homogeneous broadening (see Fig. 10). A brief 
introduction to the physics of 2D photon echo spectroscopy is given in the appendix c (for 
details see for example [19]). We want to measure Full Width at Half Maximum (FWHM) 
of the homogeneously broadened peak[22]. We consider Fig 2 (a) of G. S. Engel et. al[20]. 
The homogeneous broadening is along the main diagonal (see Fig. 10). From the scale 
given in terms of nano-meters of the figure 2(a) of G. S. Engel et. al., the FWHM is about 
~ lOnm and the exciton peak occurs at 810 nm. This gives the frequency broadening 
Sufwbm - 2-87 x 10 13 Hz. 

With this experimental information we plot I (00 — uq) such that FWHM is about ~ 
lOnm = 2.8710 13 iJ2;. Clearly for the Lorentzian, at FWHM Slofwhm = 277. This gives 
7] = 0.0143 fsec- 1 . 

C. Long coherences 

We now have all the required parameters, from the experimental information, namely, 
77 = 0.0143 fsec -1 , J = 100cm -1 , and A = 100cm -1 . With these values we plot the 
dynamics of the density matrix elements r(t), x(t), and y{t). We clearly see that the density 
matrix elements show oscillations upto 500/sec, mimicking the long coherences observed in 
the experiments of G. S. Engel et. al. To reproduce the actual spectra observed for example 
by G. S. Engel et. al. one has to go beyond this simple two state model. One has to 
consider a detailed model of the FMO complex containing not the two coupled molecules 
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ij = 0.0143 fsec - 



tj= 0.0143 fs" 




FIG. 11: Line shape function I(u) for A = 100cm ,J= 100c?n , rj = 0.0143 fsec 1 (upper 
left). This value of rj i.e., 0.0143 fsec -1 give the correct value of homogeneous line-broadening 
5uj ~ 2.87 x 10 13 Hz. 

(as considered here) but the seven coupled BChl molecules and there interactions with the 
protein matrix. This clearly requires a considerable computational challenge and one has to 
rely on the numerical approach rather than on an simple analytical solution as given here. 

IV. CONCLUSION 



We have seen that 2nd Born quantum master equation fails in modeling the real light 
harvesting systems. The reason is that it is perturbative in nature. One has to develop non- 
perturbative theories for the problem. Also resent studies with sophisticated 2D photon 
echo spectroscopy show long livid coherence effects in exciton motion which contradicts the 
long held old idea of incoherent motion [4]. New approaches should explain this by going 
beyond the perturbation theories or new mathematical framework is needed. 

Our approach towards the 2nd Born master equation enable us (1) to introduce an effi- 
cient numerical scheme, (2) to quantitatively known the validity regime of the Markovian 
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approximation. We saw that it is inadequate for the present problem. In the last section 
we introduced a phenomenological approach. In this, we modeled the effect of bath as a 
stochastic noise and its strength was calculated from motionaly narrowed lineshape. For 
this we used Kubo's stochastic theory of lineshapes. Significant role was played by the 2D 
photon echo spectroscopy as we had the important information about motionaly narrowed 
lineshape. We saw that the density matrix elements showed oscillations upto ~ 500/sec, 
thus mimicking the long coherences effects as observed in recent experiments. 
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VI. APPENDIX 

A. Novikov's Theorem 

By taking the matrix elements of equation (19) 

in ?P^ = ih ^ mm = (l\[H,p(t)}\2), (30) 

one can formally solve for pi2(t) as 

i r f 

Pn(t) - p 12 (0) = -- / ds(Ap 12 (s) 

+ J(l-2pii(s))-e(s)pi2(s)). (31) 

Taking the functional derivative ^K of the above equation wrt e(t) and plugging into 

(e(t)p l2 (t)) = h 2 V J + °° ds5(t - s) (j^-) (32) 

yields the result equation (22). 
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B. Long lime solution 



Laplace transform of the system takes the form 

( s 2J 
s + 7] -A 
V -23 A s + nl 

After inversion, in the long time limit, one has 





' r(s)^ 




( l ) 




x(s) 


= 







\y(s) j 




I - j '° ) 



t\- 



TlJ* 



r(t) ~ 1/2 + rational function(J, A, r])e v' a +j' a +^ 2 



This takes the value 1/2 when t » t re i ax 



ri 2 +J 2 +A 2 
ir,J 2 



C. Brief introduction to 2D photon echo spectroscopy 

A brief overview of 2D photon echo spectroscopy is given with an emphasis on the un- 
derlying physics of multidimensional echo spectroscopy (see for details[19]). 2D photon echo 
spectroscopy (2DPES) is a kind of generalization of the Pump-Probe Spectroscopy (PPS). 
In optical PPS, an ultra-short pump pulse with wide bandwidth creates excitation of various 
electronic transitions and the subsequent probe pulse selectively measures the transient ab- 
sorption of the electronic states. This transient probe absorption is a function of delay time 
between the pump and the probe pulses. Thus one can get dynamical information (tempo- 
ral changes of absorption) of the relaxation processes. But the pump-probe spectroscopy is 
insensitive to the coherences created by the optical excitation whereas 2DPES is coherence 
sensitive, it can temporally resolve the dynamics of the coherence (off-diagonal elements 
of the density matrix). In 2DPES three ultra short pulses are send through the sample. 
The first pulse creates the coherence state between the ground and excited statesp^e . This 
evolves for some time period r (order of femto seconds), then the next pulse interacts with 

(2) 

this already excited system. This yields either the ground state p gg or inter-exciton coher- 

(2) 

ence state p ee ,. This doubly excited state then evolves for another interval of time called 
population time T until a third pulse interacts with the system. This third interaction fi- 
nally yields 3rd order density matrix elements such as p e , which emits an echo signal (in 
phase matched direction) by decaying after a time interval t. This can be used to measure 
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real time dynamics of resonance coupling in FMO systems and this can shed light on the 
conformal changes in molecular structure such as hydrogen bound breaking (see below). 

The usual linear spectroscopic methods like linear absorption or pump-probe spectroscopy 
can only provide highly averaged information about the system under study, for example, 
in linear absorption spectra the broadening is both due to homogeneous broadening (HB) 
and inhomogeneous broadening (IHB). But 2D photon echo spectroscopy can resolve these 
two contributions. 

The way in which it resolves can be explained as follows. Consider that an ultra-short 
pulse perturbs the system at time t = 0. Consider that our system is composed of several 
chromophors with different electronic transition frequency (static in-homogeneity) and it is 
interacting with thermal bath (some protein). The first pulse creates the coherence state 
between the ground and excited state p ge of a choromophore due to the dipolar matrix 
elements coupling the ground and excited state (considering that the light-matter interaction 
is treated with first order perturbation theory-weak field regime). We have an ensemble of 
coherences p ge with a specific phase relation at t — 0. Due to static in-homogeneity the phase 
relation between these density matrix elements will be lost with time (phase randomness) 
but it will re-appear after sufficiently long time ! (if we consider for the moment that there is 
no bath and no random perturbation of the phases). Let this coherences (pge (t)) evolves for 
some time period r (order of femto seconds), then the next ultra short pulse interacts with 
this already excited and evolving system. This yields either the ground state populations 
p gg (no coherences) or inter-exciton coherence states p e J. These doubly excited states then 
evolves for another interval of time called population time T until a third pulse interacts 

(3) 

with the system. This third interaction has an opposite effect and creates the coherences p eg 
(which are complex conjugates of the first coherences). The time evolution of this exactly 
cancel the "phase randomness" developed in the initial time interval r (because evolution 
operator for p eg is the complex conjugate of the evolution operator for p ge . If there is no 
bath (system is isolated) then after an interval of time r the phases again "cohere" (they 
assume the same distribution as they had at time t = 0) and finally this "re-locking" of 
phases yield an echo signal (in phase matched direction). 

Now consider that our system is interacting with the bath (our system is open). This 
cause a "stochastic phase randomness" between the phases of p ge of various chromophors in 
the ensemble. If the population time T is sufficiently long the phase relationship between 
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density matrix elements of the chromophores will be permanently lost and and no echo will be 
seen. Thus, we can say that the maximum population time T max directly depends upon (a) 
the strength of system-bath coupling, (b) measure of the "fastness" of the bath fluctuations. 
In Kubo's stochastic theory these are respectively A and 7. Thus T max tell us about the 
character of homogeneous broadening mechanism. The t — t correlation or in frequency 
domain uo T — u T correlation (along the diagonal direction in 2D spectra) tell us about the 
in-homogeneous broadening, as the experiments are done by systematically varying r and 
at a given r (the time gap between the first two pulses) the amplitude in the 2 D spectrum 
after that given time r from the second pulse (i.e., after population time) show a correlation 
in the form of the elongation of the peak in the diagonal direction, a direct signature of 
in- homogeneous broadening. Thus one get the information about both HB and IHB. 

2D photon echo spectroscopy can also tell us about the real time resonance coupling dy- 
namics and bond breaking dynamics (Fig. 12). Consider that we have two tuning forks with 
characteristic frequencies ui\ and u>2- Let us excite this system with a spectrum of frequen- 
cies and detect the amplitude of vibration with some frequency analyzer (some electronic 
instrument) and plot various frequencies along the x-axis (exciting frequencies) and y-axis 
(detecting frequencies). We will see two peaks occurring at Ui and 101 in the "2D spectrum" 
(along the diagonal Fig. 12). 

Now consider that our tuning forks are coupled by some spring (say we have weak cou- 
pling). Then again repeat the experiment and plot the 2D spectrum. This time we will see, 
along with the diagonal peaks, two "cross-peaks" along the anti-diagonal direction. These 
cross peaks are the consequence of coupling. If we further reduce the spring coupling the 
magnitude of these cross peaks diminish and finally disappear with our removal of coupling 
springs. 

This mechanical analogy can be directly applied to the changes in the molecular structure. 
In Fig. 1 hydrogen bound breaking dynamics between two chemical groups is shown. Two 
chemical groups have characteristic frequencies oj\ and U2- Two peaks appear at these 
frequencies in the 2D photon echo spectrum. Two cross-peaks also appear as shown which 
is due to the coupling of these two chemical groups. The coupling is due to the hydrogen 
bound. When the bound is broken this cross-peaks also disappear. This can be tracked 
in real time by varying the magnitude of the population time. Similarly in the electronic 
spectrum the cross peaks tell us about the electronic coupling between the chromophores 
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FIG. 12: (A) 2D spectroscopy can track (in real time) the changes in the molecular structure. In 
the above example hydrogen bound breaking dynamics between two chemical groups is shown. Two 
chemical groups have characteristic frequencies io\ and W2- Two peaks appear at these frequencies 
in the 2D photon echo spectrum. Two cross-peaks also appear as shown which is due to the 
coupling of these two chemical groups. The coupling is due to the hydrogen bound. When the 
bound is broken this cross-peaks also disappear. (B) A mechanical analogy: The occurrence of 
cross peaks in the exciting-response frequency spectrum is due to the spring coupling between two 
tuning forks of frequency wi and u)2- The amplitude of cross-peaks in the 2D spectrum measures 
the strength of coupling. When the spring constant decreases the cross peaks diminishes and finally 
goes to zero with no coupling. 

(the system-system coupling, usually denoted as J or V) also called resonance coupling. 
The oscillation in the magnitude of the cross peaks show the oscillations of the pigments 
about their equilibrium positions (in physical space) as the variation about the equilibrium 
position modulate the resonance coupling strength. 

A brief mathematical formulation can be described as follows. In semi-classical approxi- 
mation for field-matter interaction the interaction Hamiltonian is written as 

H int = -/2.E(r,t) (33) 

Here dipole approximation (weak variation of the electric field amplitude over the size 
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of the pigment) is used. If the magnitude of the above Hamiltonian is much weaker than 
the magnitude of the pigment Hamiltonian H P i g — H e + H p h + H e _ p h, then the Hi nt can be 
treated as a perturbation. 

The dynamics of the total system (the pigment) can be described by Liouville-von Neu- 
mann equation for the density matrix p. 

^ t =- l j\H ptg + H mtl p] (34) 

Now consider that the sample is interrogated with three consecutive (in time) ultra short 
laser pulses. Treating H int as a perturbation the third order density matrix is given as 

/*oo /*oo /*oo 

p {3) (r,t) = (-) 3 dt 3 dh dti6(t 3 )e-^x (35) 

" Jo Jo Jo 

^e(t 2 )e-*° 2 L M e(t 1 )e-^ 0l L M p(-oo)E(r, t - t 3 )E(r, t - t 3 - £ 2 )E(r, t-t 3 -t 2 -t x ). 

Here C* = [H pig ,*] and L p * = [//,*]. Experimentally we do not measure pit) but the 
induced polarization p( 3 \v,t) = tr{pp( 3 \t)}. 

poo poo /*oo 

P (3) (r,t) = (-) 3 / dt z dt 2 rft 1J R^(t 3 ,t 2 ,ti)E(r,t-t 3 )E(r,t-t3-t 2 )E(r,t-t3-t 2 -t i ; 

(36) 
The 3rd order response function is given as 

R i3) (t 3 ,t 2 ,h) = (^) 3 e(t 3 )e(t 2 )e(t 1 )([[[/i(t3 + t 2 + t 1 ), / i(t 2 + t 1 )], / i(t 1 )], / i(o)p(o)]). (37) 

The time dependence of /i can be shifted to the time dependence of the density matrix 
using time evolution operators (resulting p, as time independent). Ultimately one have to 
take the trace over the bath modes thus one will end up having an expression involving time 
depended reduced density matrix elements of the system. 

The electric field of the incoming pulses can be given as 



E(r, t) = Y^ (tiiEiit)^^ 1 *-^ + c.c.) (38) 



i=l 



Where n is the unit vector pointing in the direction of the field and E(t) is the temporal 
envelope of the pulse. The non-linear polarization can also be expanded in the Fourier 
components 
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p( 3 ) ( r , t) = Y^ Pj 3) (t)e ik - r "^ 4 (39) 

i 
Where 

k ; = ±ki ± k 2 ± k 3 , u>i = ±Wi ± w 2 ± w 3 (40) 

This non-linear polarization will emit a signal electric field in the phase matching direc- 
tion. It can be shown?? ( if the electric field amplitude varies slowly over the pigment size 
and if the refractive index of the material is frequency independent) that 

^nal OC P (3) (41) 

Now under the impulsive limit (excitation femtosecond laser pulses as delta functions) 
Ei{t) = E x 8(t + T + r), E 2 (t) = E 2 8(t + T), E 3 (t) = E 3 S(t) the signal electric field 
becomes proportional to the response function. The Double Fourier transform of the signal 
field gives the complex 2D spectrum 

/+oo /"+oo 

dt dr^(t,T,r)e^ e ^ (42) 

-oo J — oo 

This is the 2D photon echo signal. The whole problem boils down to the calculation of 
3rd order response function which further involves the time evolution of the reduced density 
matrix. 
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